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EXPLICIT-IMPLICIT SCHEME FOR RELATIVISTIC RADIATION HYDRODYNAMICS 

HiROYUKi R. Takahashi\ Ken Ohsuga^, Yuichiro Sekiguchi Tsuyoshi Inoue^, and Kengo Tomida^' ^ 

ABSTRACT 

We propose an explicit-implicit scheme for numerically solving Special Relativistic Radiation Hy- 
drodynamic (RRHD) equations, which ensures a conservation of total energy and momentum (matter 
and radiation). In our scheme, 0th and 1st moment equations of the radiation transfer equation are 
numerically solved without employing a flux-limited diffusion (FLD) approximation. For an hyper- 
bolic term, of which the time scale is the light crossing time when the flow velocity is comparable to 
the speed of light, is explicitly solved using an approximate Riemann solver. Source terms describing 
an exchange of energy and momentum between the matter and the radiation via the gas-radiation 
interaction are implicitly integrated using an iteration method. The implicit scheme allows us to relax 
the Courant-Friedrichs-Lewy condition in optically thick media, where heating/cooling and scattering 
timescales could be much shorter than the dynamical timescale. We show that our numerical code 
can pass test problems of one- and two-dimensional radiation energy transport, and one-dimensional 
radiation hydrodynamics. Our newly developed scheme could be useful for a number of relativistic 
astrophysical problems. We also discuss how to extend our explicit-implicit scheme to the relativistic 
radiation magnctohydrodynamics. 

Subject headings: hydrodynamics - radiative transfer - relativity 



1. INTRODUCTION 

Relativistic flows appear in many high-energy astro- 
physical phenomena, such as jets from microquasars 
and active galactic nuclei (AGNs), pulsar winds, mag- 
netar flares, core collapse supernovae, and gamma-ray 
bursts (GRBs). In many of these systems, the mag- 
netic field has a crucial role in dynamics. For exam- 
ple, magnetic fields connecting an accretion disk with a 
central star, or, different points of accretion disks are 
twisted and amplified due to the differential rotation, 
launching the jets ( Lovelace 19761: IBlandford fc Pavnd 
[TOSllUchida fc Shibata 1985). The production of astro- 
nomical jets has been observe d in ma g netohydrodynamic 
(MHD) simulations (iHavashi et~all Il996t iMeier et all 
120011: iKato et all |200j). The jets are also powered by 
a subtraction of the rotation al energy of a central black 
hole via th e magnetic field (Blandford & ZnaieW 119771 : 
iKoide et a l. 2002; McKinncy & Blandford 2009). Also 
the magnetic fields play an important role in accre- 
tion disks to transport the angu lar momentum outward , 
leading to the mas s accretion (jBalbus fc HawlevI 119911 : 
IHawlev et all[l995h 

Not merely the observational point of view, but the 
radiation field is also an important ingredient in the 
dynamics of relativistic phenomena. The radiation 
press ure force would play a key role in jet accelera- 
tion ( Bisnovatvi-Kogan fc Blinnikovl Il97"7t iSikora et al.l 



[1991: lOhsuea et all 120051 : lOkuda et al. llooi). Recently 
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iTakeuchi et al.l (j2010f ) showed a formation of radia- 
tively accelerated and magnetically coUimated jets us- 
ing non- relativistic radiation MHD (RMHD) simulation s 
(see also lOhsuga et al.ll2009l : lOhsuga fc Mineshigell20lH) . 
The radiative effect also plays important roles in non- 
relativistic phenomena. For example, the dominance 
of the radiation pressure in optically thick, geometri- 
cally thin (standard) accretion disk s could be unstable to 
thermal a nd viscou s instabilities (Lightman fc Eardle 
1974: Shibazaki fc Hoshi 19 75: Sliakura fc Sunyacv 197j' 



Takahashi fc Masadall201ll ). Here we note that the global 



radiation hydrodynami c (RHD) simulation s showed the 
growth of instabilities (jOhsugal [2006) , but IHirose et al.l 
pOO&) concluded that the radiation pressure dominated 
region is thermally stable using RMHD simulations of a 
local patch of the disk. Thus the high energy phenomena 
would be understood only after we take into account the 
magnetic field and the radiation consistently. 

Many methods have been proposed to include radiative 
effects in non-relativistic plasmas. Solving a full radia- 
tion transfer equation is the most challenging task for di- 
rectly connecting results between nu merical simulations 
and observations (jStone et al.lll992f ). But it is a hard 
task to couple with the hydrodynamical code due to its 
complexity and high computational costs. Another class 
of simple method is based on the moment method. In the 
Flux-Limited Diffusion approximation (FLD), a 0th mo- 
ment equation of the radiation transfer equation is solved 
to evaluate the radiation energy density, E'^,. The radia- 
tive fiux, F'j. , is given based on the gradient of the radia- 
tion energy density without solving the 1st moment equa- 
tion. This is quite a useful technique and gives appro- 
priate radiation fields within the optically thick regime, 
but we should keep in mind that it does not always give 
precise radiation fields in the regi me where the optical 
depth is around unity or less (see lOhsuga fc Mineshigd 
I2011f ). Thus it is better to solve the 0th and 1st moment 
equations. 
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When we solve two moment equations, the closure 
relation, F/*^ — D^^E'^,, is used to evaluate the radi- 
ation stress, P^'-', where D^^ is the Eddington tensor. 
In the Eddington approximation, the tensor is simply 
given by D'-' = 6^^ /3, where 6^^ is the Kronecker delta. 
This implies that the radiation field is isotropic in the 
comoving frame, which is a reasonable assumption in 
the optically thick media. Another class of the Edding- 
ton tensor is obtained by approximately taking into ac- 
count the anisotropy of radiation fields (so called M-1 
closure, .Kershay 19^ 6; Mincrbo 1978; Lcvcrmorc 1984) • 
iGonzalez et al.l ()20d7f ^ implemented the M-1 closure in 
their non-relativistic hydrodynamic code and demon- 
strated that the anisotropic radiation fields around an 
obstacle are accurately solved. 

In the frame work of relativistic r adiation hydrodynam- 
ics (RRMHD), iFarris et af] (|2008[ ) first proposed a nu- 
merical scheme for solving general relativistic radiation 
magnetohydrodynamic equations a ssuming an optically 
thick medium (i.e., D'^ = 5'^ /i). iZanotti et a"IT(|201l( ) 
implemented the general relativistic radiation hydrody- 
namics code and adopted it for the B ondi-Hoy le accretion 
on to the black hole. Recentlv. Shibata et al.l ([2011,) pro- 
posed another truncated moment formalism of radiation 
fields in optically thick and thin limits. 

In the RR(M)HD, there are three independent 
timescales. The first one is the dynamical time scale 
tdyn = L/v, where L and v are the typical size of the 
system and the fluid velocity. The second one is the 
timescale at which the characteristic wave passes the sys- 
tem, tu, = L/\, where A is the characteristic wave speed. 
These two timescales are comparable and close to the 
light crossing time, L/c (where c is the speed of light) in 
relativistic phenomena (i.e., w ~ A ~ c). The third one 
is the heating/cooling and scattering timescales, tab and 
tsc- These timescales become much shorter than the dy- 
namical one in the optically thick medium, since the gas 
and the radiation are frequently interact. This indicates 
that equations for the radiation fields become stiff. In 
the explicit scheme, a numerical time step Ai should be 
restricted to being shorter than these typical timescales 
(At < min[trfy„, tu,, tab, tsc]) to ensure the numerical sta- 
bility. Therefore, high computational costs prevent us 
from studying long term evolutions when the gas is opti- 
cally thick, if we are interested in the phenomena of the 
dynamical timescale. 

In this paper, we propose an explicit-implicit numerical 
scheme overcome this problem. As the first step, we con- 
sider the relativistic radiation hydrodynamics by neglect- 
ing the magnetic field (but, see section [5] for the discus- 
sion how to extend our numerical scheme to RRMHD). 
The proposed scheme ensures a conservation of total en- 
ergy and momentum. Governed equations are integrated 
in time using both explicit and implicit schemes. The for- 
mer solves an hyperbolic term (timescale ^ tdyn)- 
The latter treats the exchange of energy and momen- 
tum between the radiation and the matter via the gas- 
radiation interaction, whose time scales are ~ tab and 
tsc- This method allows us to take a larger time step 
At > tab, tsc than that of the explicit scheme. 

This paper is organized as follows: In § [5J we intro- 
duce argument equations for RRHD and the numerical 
scheme is shown in § [3l Numerical results of one- and 



two-dimensional tests are shown in § 2] Discussion and 
summary are appeared in § [5] and § [6l 

2. BASIC EQUATIONS 

In the following, we take the light speed as unity. 
The special relativistic radiation magnetohydrodynamic 
equations of ideal gas consist of conservation of mass, 

{pu").. = 0, (1) 
conservation of energy-momentum, 

(I^^D+T^rd),. = 0, (2) 

and equations of radiation energy-momentum 

TZ.. = -G^ (3) 

where p is the proper mass density, — 7(1, w') is the 
fluid four velocity, and T^J^ and T^^^ are energy momen- 
tum tensors of fluid and radiation. Here 7 = vT+~Ui?F 
is the bulk Lorentz factor and is the fluid three veloc- 
ity. Greek indices range over 0, 1,2,3 and Latin ranges 
over 1,2,3, where indicates the time component and 
1, 2, 3 do space components. 
The energy momentum tensor of fluids is written as 

T^^^p^u^u'^+Pyil^", (4) 

where Pg is the gas pressure and diag 77^" = (—1,1,1,1) 
is the Minkowski metric. The specific enthalpy of rela- 
tivistic ideal gas 1^ is given by 



where F is the specific heat ratio. 
The energy momentum tensor of radiation is written 

as 

^ = (^;^) 

where Er, F^, and P^^ are the radiation energy density, 
flux and stress measured in the laboratory frame. 

The radiation exchanges its energy and momentum 
with fluids through absorption/emission and scattering 
processes. The radiation four force is explicitly given 

by 

G" = -pK (47rB7 - jEr + UiF^.) 

~ pas [WEr + lu^UjP^^ - (72 + «2) , (7) 

and 

= -AttpkBu' + p{k + crs)(7F; - u^P;^') 

- posu' (72^, - 2-iUjF} + UjUkPt) . (8) 

where u — y/uiU"^. k and CTs are absorption and 
scattering coefficients measured in the comoving frame. 
Thus, equation Q is a mixed-frame radiation energy- 
momentum equation that the radiation field is defined in 
the observer frame, while the absorption and scattering 
coefhcients are defined in the comoving frame. 

In equations ([7]) and ([S]), we assume the Kirchhoff- 
Planck relation so that the emissivity ry is replaced by 
the blackbody intensity B as ry = kB in the comoving 
frame. 



Relativistic Radiation HD 



3 



The blackbody intensity B is a function of gas temper- 
ature T as B = aflr"^/(47r), where is related to the 
Stefan-Boltzmann constant crgB = a_R/4. The plasma 
temperature is determined by 



Pa 



(9) 



where fc^ and nip are the Boltzmann constant and the 
proton mass, and is the mean molecular weight. 

Since we consider the moment equation of radiation 
fields, we need another relation between Er, and P^^ 
to close the system, i.e., the closure relation. Here we 
leave it as a general form described by 



(10) 



where dash denotes the quantity defined in the comoving 
frame. The explicit form of closure relation is introduced 
in § 13.31 after the formulation. 

In the following, we deal with the one-dimensional con- 
servation law in the x-direction: 



'dt 



dx 



(11) 



which follows equations (HI, (HI, and ([3]). Primitive vari- 
ables V are defined as 



V = 



u 
Pa 



(12) 



Conserved variables fiuxes J^, and source terms S have 
forms of 



U 



Et 



(13) 



r 

pxk 



xk 



(14) 



and 



S 



1 ° ^ 




( ° \ 


















Se 




-G" 









(15) 



where _D = p7 is the mass density measured in the labo- 
ratory frame. The total energy density Et and momen- 
tum density mt are given by 



Et = i^HD + Er = pCT^ -Pg+Er 



and 



F^ 



(16) 



(17) 



where i?HD and to^j-, denote the energy and momentum 
density of the fluids, respectively. 

It should be noted that Er and Fr are not only the 
primitive variables, but also the conserved variables. 



Thus, it is straightforward for the radiation fields to con- 
vert the primitive variables from the conserved variables, 
as we will see in S 13.21 



3. NUMERICAL SCHEME FOR RRHD 

In this section, we propose a new numerical scheme for 
solving RRHD equations. A conservative discretization 
of 1-dimensional equation (|lip over a time step At from 
t ~ n/St io t ~ (n + l)At with grid spacing Ace is written 

as 



n+1 



At 
Aa- 



(/.. 



(18) 



where is the conservative variable at a; = Xi and t — 
nAt, and is the upwind numerical fiux at the cell 
surfaces x = x^j.i. In the numerical procedure, we divide 
equation (|18p into two parts, the hyperbolic term 



and the source term 



(19) 



(20) 

where lA* is the conservative variable at the auxiliary 
step. In the following subsection, we describe how to 
solve these two equations. 

3.1. hyperbolic term 

For the hyperbolic term, equation is solved as an 
initial value problem with the initial condition 



Z^(r,a;) 



Ul ior X < x^^i 
Uji for X > .Tj^i 



(21) 



where the subscript L and R denote left and right con- 
stant states on the cell interface. When we take Ul — l^i 
and IAr = tii+i^ the numerical solution is reliable with a 
first order accuracy in space. 

In order to improve the spatial accuracy, primitive vari- 
ables at the zone surface Vi±i/2 are calculated by inter- 
polating them from a cell center to a cell surface (the 
so-called reconstruction step). The primitive variables 
at the left and right states are written as 



2 ' 



(22) 



where we take S = L{R) with the plus (minus) sign. 
The slope 5xP should be determined so as to preserve 
the monotonicity. In our nu merical c o de, w e adopt the 
harmonic mean proposed by Ivan Leeil (|1977f) , which has 
a 2nd order accuracy in space, as 



2max(0, AP+AT^-) 
AV+ + A-P_ 



(23) 



where 



AV± = ±{V.,±i - V^) 



(24) 

Here we adopted a 2nd order accurate scheme, but 
a highe r order scheme such as Piecewise Constant 
Method (E omissaro ^[19991) Piecewise Parabolic Method 
(|Colella fc WoodwardI 11981 IMarti' fc MueileJ [19961) or 
the other scheme can be applicable to both hydropart 
and radiation part. 
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By computing the primitive variables at the cell sur- 
face, the flux ^i±i/2,5 is calculated directly from 'P^j-i g 

(but see § 13.31 for the procedure of calculating P^'^ ap- 
peared in equation [T3] from P^^^). Then, numerical fluxes 
are calculated using an app roximate Riemann solver. We 
adopted the HLL method (Hart en et al.lll983[) . which can 
capture the propagation of fastest waves. The numerical 
flux is then computed as 



equations ([28)) - ([29]) become non-linear equations for 
jjn+i ^ Another difficulty comes from the closure relation 
ptj = D'^Er. The Eddington tensor D^' = D'^{Er,F^) 
is generally a non-linear function of Er and F^, so we 
cannot adopt the simple implicit method by replacing 
Er El'+^ and F^ F^+^'' in the Eddington tensor. 

In this paper, we propose an iterative method, which 
solves the following equation 



A+J-,+ i,i - A_J-,+ i,^ + X+X-{U,+ i^,n-K+lL) Ei"^+^^ = e; + ^tSE{Ei^+'\F^^+^^'\P^^+^^'=\V^^ 



(r, 



A+ - 



(25) F 



Here A_ and A+ are 

A_ = min(0,AL-, A_R_), 



(26) 



and 



A+ = max(0, \l+,\r+), 



(27) 

where A5+ and As_ are the right and left going wave 
speed of the fastest mode. For example, they corre- 
spond to the sound speeds for the relativistic hydro- 
dynamics. When the radiation field is included, the 
light mode determines the fastest wave speed, which de- 
pends on the closure relation. The fastest wave speed is 
discussed after specifying the closure relation in § 13.31 
We note that although we adopted the HLL scheme 
for simplicity, the higher order approxi mate Riemann 
solvers such as the HLLC (iMignone fc Bo do 2005, 2006; 
Honkkila fc JanhunenI 120071 ) and HELD (|Mignone et al.i 



2001 " can be implemente d in the relativistic ra diation 



(Magneto)hydrodynamics (ISekora fc Stonel[2009D . 
Using the numerical flux /_(_ 1 , the conserved variables 

U* is obtained from equation (|19p . It should be noted 



that the total energy density E^'^^ and momentum den- 
sity m"^^' ~dX t = {n-\- l)At are already solved although 
we only integrate the hyperbolic term. Thus, when the 
radiation energy density and flux are obtained by solv- 
ing equation (1201) . the fluid energy density and momen- 
tum density are immediately computed (see, equations 

mMi). 

3.2. Source term 

Next, we show how to solve equation ([20)1 . The source 
term appears in equation ([3]) , which treats the interaction 
between the radiation and the matter. As discussed in 
§[1] the heating/cooling and scattering timescales can be 
shorter than the dynamical timescale in optically thick 
media. It prevents us from studying the long time evolu- 
tion (^^ tdyn) when the equation is explicitly integrated. 
To overcome this difficulty, we want to construct an im- 
plicit scheme as 

E';+^ = e; + A^5i^(£;;+^F,"+l^^P,"+l■^■^P,^+l), 

(28) 

^ f:^' + MS'p{E';+\Fi'+^^' ,p^+^'^\v]i+^), 

(29) 

where Vh is primitive variables of fluids (i.e., p, u^, pq) . 
We confront two problems to solve equations ([?5|) - (^9)) . 
The first problem comes from the appearance of VJ^^^- 
Since PJ^~^^ should be obtained after computing 



(m+l),i p'^^i 



where m = 0, 1, 2... indicates the iteration step. We take 



E, 



(0) 



El\ F 



(0)J 



F"'^, and P 



(0) 



PJI^ for the initial 



guess (also E* and F*'^ can be candidates for the ini- 
tial guess). We note that the hydrodynamic quantity is 
explicitly evaluated at m-th step due to the complexity 
discussed above. 
Next we introduce the following two variables 



(32) 
(33) 

By substituting equations p2p and p3p into equations 
([50]) - ([5T|). and taking the first order Taylor series in 

(5£;^™+^^ and we obtain 



5E., 



5Fr 



(m+l)j 



^'"+^^ ^ / Ats'iP^ + E* ~ E^r^ \ . 

, \Aisf^'' -I- f;^^ - p'r^'^f'^ ' 

where S^^^ = SE[Ei"'\ F^"'^'', P^^^^^'^P^'^^] and 
4")^* = 5^[£;^™^i^,('")'■'■,P,(™)^^■^pM]. Here C is the 
4x4 matrix given by 



C 



where 




dSE 
dEr 
OSe 
dP} 

dP^i 
dSi 



dP''' asE bse I dP''' asE 

dEr dP''' ' dF^ dF^ dP''' 

dP;±dS}_ dsj dP''' dS'p 

OEr dP''' ' OF' dF^ OP''' 



---npj + asp-fu , 
-asPJUkUi, 



dEr 

dP} 
dS^F 
dPr^^ 



^=<JsPl'u\ 



-Kp'-fS^j — aspj {Sj + 2u^Uj) 



(35) 

(36) 
(37) 
(38) 
(39) 
(40) 



= Kp5].Ul + (Tsp [Slui + u'-UkUl) . (41) 



Here we use equations ([7]), ([S]), and (ITSI) . Also 
{dPr^)/{dEr) and {dPr')/{dFr ) are required to com- 
plete matrix elements. These quantities depend on the 
closure relation given in equation (|T0| . We do not specify 
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the closure relation here, but the explicit form of these 
quantities is shown in § 13.31 and appendix El 

When P^^ is the linear function of Er and F* (e.g., the 
Eddington approximation), equations pO|) - reduce 
to 



E. 



(rn+l) 



F, 



{m+l)A 



El^ 



s^^\e, 



= o,f; = o)Ai 
(£;,. = o,f; = o)At 



(42) 

By inverting the 4x4 matrix C in equation (IMl) or (|^^ . 
we obtain the radiation energy and flux at (m + l)-th it- 
eration step from equations p2)) and (|33| . In general, 
the matrix inversion is time-consuming and sometimes 
unstable. The matrix C is however, only 4x4 matrix 
so that we can invert it analytically. We also tried in- 
verting it using LU-decomposition method. We obtain 
the inverse matric stably in both scheme, but the 
analytical method is faster than the LU-decomposition 
method. Thus we decide to use analytical expression of 

c-\ 

Next, we calculate the primitive variables 'p(™+^) from 
updated conservative variables U'^"^~^^\ As pointed in 
§ [21 i?r and F!^ are the both conserved and primitive 
variables. Thus, the recovery step is unnecessary for 
the radiation field. We need to compute 'P^™^^-' from 

^(m+l)_ 

Since the total energy E^^^ and the momentum 
^n+i,/c already determined, the energy and momen- 
tum of fluids (i?HD^^'' "^HD^^^ '^) can be calculated 



(m+l) 



(m+l),fe 
^HD 



n+l,k 



^(m+l),fc 



(43) 
(44) 

Then, three unknown variables p^'^^^\u^™-^^^-^ .■p'g^^'^^ 

are computed from D", m[j"^+^)''-, and 4™+''. Thus, 
the recovery step in RRHD is the same as that in rela- 
tivistic HP. W e adopt the recovery method developed by 
iZenitani et al.l (2009) for solving a quartic equation. We 
briefly show the method. In the following discussion, we 
drop the superscripts n and (m -|- 1) for simplicity. 

The gas density _D, momentum to^q, and energy E'hd 
are related to the primitive variables p, as 

B = p7, (45) 
{P + Tipg)ju\ (46) 

^{p + r,Pg)^', (47) 

where Fi = r/(r — 1). From equation P^l) . we obtain 



"■HD = 



Pa 



1 / |?^hd| 

7ri I 



D 



(48) 



where and rriHD = Y^u ^^hd^hd- By substituting equa- 
tion (j48l) into (|47l) . we obtain a quartic equation on 

U = \JUiV} 

J{u) = Tl (Sgo - ml^) - 2rimHDi?M' 
+ [TlElj, -D^- 2ri(ri - l)|mHDn 



-2(ri-i)i?|mHD|w-(ri 



0. (49) 



By solving above quartic equation, pg, p and are 
computed from equation (|48l) . p5l) . and p6t . When 
we solve eq uation (1491) . we first adopt a Brown method 
(jNunohiro &: Hirano 2003), which gives analytical so- 
lutions of a quartic equation. If reasonable solutions 
are not obtained, we solve equation (j^^ using Newton- 
Raphson method with an accuracy of f{u) < 10~®. It is, 
however, noted that numerical solution converges with- 
out switching to Newton-Raphson meth od in the test 
problems shown in § |4l Also we note that IZenitani et al.l 
(,2009) showed that the Brown method can be applicable 
to obtain solutions with a wide range of parameters. 

When all the primitive variables are recovered, we ob- 
tain solutions pt'^+i). These solutions would satisfy 
equation (|42l) . But we note that while E^ , and P^^ 
are evaluated at (m -I- l)-th time step, we evaluate Vh 
at (m)-th time step when we solve equation The 
evaluation of Vh at (m)-th iteration step might be prob- 
lematic when the density or temperature jump is large 
(e.g., at the shock front). Thus, we again solve equa- 
tion (|42l) using updated primitive variables vlm + 1) and 
check the difference between two successive variables, 
\6E^"'+^^\/Ei"'+^\ 



^^p(m+l),^^^ and 



(m) I 



,(m+l)| 



When these quantities is larger than a 



specified value (typically ~ 10~^), the source term is in- 
tegrated again using updated primitive variables. This 
process is continued until t wo successive variables fall be- 
low a specified tolerance (jPalenzuela et al.1l2009[ ). Sim- 
ilar iteration method is also adopte d in relativis tic re- 
sistive magnetohydrodynamics (Pale nzuela et all [2009). 
Since Vh is evaluated at the current iteration step, so- 
lutions obtained by solving equation (|42p might not be 
converged when the shock appears. But as we will see in 
§ m we obtain solutions correctly even when the cooling 
time is much shorter than At (see, $ I4.1.1|) . or a shock 
appears (see, § l4.2p . 

Now we summarize our newly developed scheme for 
solving RRHD equations. 

1. Calculate primitive variables at the cell surface 
^^±1/2,8 from Vi. 

2. Compute the numerical flux fi±i/2 using approxi- 
mate Riemann solvers (e.g., HLL scheme). 

3. Integrate the hyperbolic term from numerical flux 
/i±i/2 (in equation [TOl) . Then we obtain the inter- 
mediate states of conservative variables U* . 

4. Compute matrix elements of C given in equation 
(f35ll. Then, calculate F^"+^^ and F^"'^^'''" by in- 



verting the 4x4 matrix. 
5. Calculate primitive variables P^^+i) from E. 



m^^+^'>''' and D"+\ 



(m+l) 
HD ' 



6. When the difference between two successive val- 
ues does not fall below a specified tolerance, re- 
peat from step 4. The updated primitive variables 



of fiuids V 



(m+l) 



are used to evaluate matrix C. 
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3.3. Closure Relation 

In the above discussion, we do not specify the closure 
relation, but take the general form given in equation 
(jlOp . When we compute the numerical flux using the 
HLL scheme, the wave speeds of the fastest mode A± are 
needed, which depend on the closure relation. Another 
term depending on the closure relation is P^^ , which ap- 
pears in numerical flux (ITil) . Also we need to evaluate 
its derivatives, dP^^ /dE,. and dP^^ /dF^ to compute the 
matrix C. In this subsection, we show how to obtain 
these quantities by specifying the closure relation. 

Many kinds of closure relation are proposed by au- 
thors as discussed in § [T] As the first step of developing 
the RRHD code, we hereafter restrict our discussion to 
the closure relation being provided when we assume the 
Eddington approximation. Then, the closure relation is 
described by 

p;^^ = — i?;, (50) 

which is valid when the radiation is well coupled with the 
matter so that the radiation field is isotropic in the co- 
moving frame. This corresponds to the Eddington factor 
being 1/3. 

First, we show the wave speed of the fastest mode. By 
assuming the relation given in equation (|50p . the char- 
acteristic wave velocity of fastest mode in the comoving 
frame is 1/ V3, which is equivalent to the sound speed in 
the relativistic regime. Thus the maximum wave velocity 
is always l/\/3 in the radiation hydrodynamics with the 
Eddington approximation. The wave velocity in the ob- 
server frame A_ and A+ is obtained by boosting ±l/\/3 
with the fluid velocity. 

Next, we show how to obtain radiation pressure P^^ 
measured in the observer frame. By performing Lorentz 
transformation on the radiation energy momentum ten- 
sor and combining it with equation (j50l) . the radiation 
stress tensor in the observer frame obeys the following 
equation 
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UkUmPr 



where 



1 + 7 



(1 + 7)^ 



(51) 



1 + 7 



u'w>UkF^, (52) 



(e.g., iKato et all [20081) . Since the radiation stress P^^ 
is a 3 X 3 symmetric matrix, we need to solve 6th order 
linear equations to compute P^^ as 

A{u)p = r, (53) 
where = iP^\ P^\ P^^ P^\ P^^ P^^) , = 

(i?", i?22^ i?33^ pl2^ ^13^ p23) g^^^ ^ ^ ^(-^) jg the 6 X 6 

matrix. Since A is a function of velocity, the derivatives 
of P^^ are described by 

. , . dp dr . , . dp dr , , 
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Fig. 1. — Thermal evolution of radiation energy E^. Crosses 
and squares denote results of explicit and implicit schemes, while 
solid curves do analytical solutions. The filled circles with dotted 
curves show the results of the implicit scheme with the time step 
of At = lOtab- 

Thus, P;^ dPi.' /dEr and dP^^/dF^ are computed by 
inverting matrix A. Since A is a 6 x 6 matrix, it is diffi- 
cult to obtain using analytical expression. Thus, we 
use LU-decomposition method to invert matrix A. Ex- 
plicit forms of A, (dr)/ (dEr), and {dr)/{dF^) are shown 
in appendix [XI 

4. TEST PROBLEMS 

In this section, we show results of some numerical tests 
for our RRHD code. This section consists of numerical 
tests for the radiation field (§ 14. 1|) and relativistic radia- 
tion hydrodynamics (§ 14. 2p . We assume that the closure 
relation is given by equation ([50]) so that the wave speed 

of radiation field is 1/ VS. 

4.1. Numerical tests for Radiation Field 

In this section, we show results of numerical tests for 
solving radiation fields given in equations ([3|). We assume 
that the fluid is static and uniform for simplicity. We 
recover the light speed as c in this subsection. Equations 
are solved with 2nd order accuracy in space. 

4.1.1. Radiative heating and cooling 

This test has been proposed bv lTurner fc Stoiii (|2001f) . 
We evaluate the validity for the integration of the source 
term appeared in equation (j20p . which is implicitly and 
iteratively integrated in our numerical scheme. For this 
purpose, we start from a static and 1 zone fluid (i.e., 
number of grid points — 1), which is initially not in 
the thermal equilibrium with the radiation. From these 
assumptions, the radiation field obeys the following equa- 
tion 

dE 

—!-= pK (47rB - cEr) , (55) 
dt 

which is analytically integrated by assuming p, k and B 
are constant, and we obtain 



47r (An 
Er{t) = — B - ( — B - Po I e 



■pKCt 



(56) 



where Eq is the radiation energy density att — 0. The ra- 
diation field approaches that in the Local Thermal Equi- 
librium (LTE, Er = 47rB/c). 
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Fig. 2. — Time evolution of Er (thick solid) and Fr (solid). The 
optical depth of the system size is 0.01. Dashed curves denote the 
light head = ct/\/3 at t = 0.5, 1.0, 1.5tc, where = L/c. Dot- 
ted curve s sh ow analytical solutions assuming steady state given in 
equation II57I I. 
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Fig. 3. — Time evolution of Er (thick solid) and Fr (solid). The 
optical depth of the system size is unity. Dashed curves denote 
the light head Ic = ctjs/Z at 4 = 0.5, 1.0, 1.5tc, where Tc = Ljc. 
Dotted curv es s how analytical solutions assuming steady state given 
in equation l|57|l . 
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Fig. 4. — Snapshot of Er at t = 1.5tc. The solid, dashed and 
dotted curves denote contours at 10%, 0.1%, 10"^% of its maximum. 

The mass density is set to p = 0.025 g cm~'^ and the 
opacity is k = 0.04 cm^ g~^, so that the corresponding 
absorption time scale is tab = ^/{pi^c) = 3.3 x 10^^ s. 
We examine two models of the initial radiation energy 
density, = IO^-Elte and 10~^i?LTE, where E'lte = 
anT^ = 10^° erg cm"'^. The specific heat ratio and the 
mean molecular weight are F = 5/3 and /.t = 1.0. 

Figure [1] shows the time evolution of Er- Squares and 
crosses respectively denote results with explicit and im- 
plicit scheme for integrating the source term with the 
time step of At = O.l^ab, while solid curves do an- 
alytical solutions obtained from equation ([55)1 . The 
data point for the model with At — O.ltab is reduced 
for plotting. The lower and upper plots correspond 
to Er{t = 0) = 10~^£'lte and lO^i^LTE, respectively. 
We can see that results with both schemes excellently 
agree with analytical solutions when the time step is 




Fig. 5. — Spatial profiles of Er ai t = 1.5rc. Asterisks, squares, 
and filled circles indicate results along y = 0, x = 0, and y = x, 
respectively. 



smaller than the absorption time scale At < tab- When 
At — lOtab, solutions with implicit scheme (filled cir- 
cle) stably reach the thermal equilibrium state, although 
they deviate from the analytical solutions before being in 
LTE. We note that the explicit scheme does not converge 
to the analytical solution when we take such a large time 
step. We also performed simulations with a much larger 
time step At = lOHab and found that solutions con- 
verge to analytical solution (not shown in figure). We 
also note that a number of iteration step used in the im- 
plicit scheme is less than or equal 2, even if At = lO^tab- 
Thus, the implicit scheme has a great advantage when 
the cooling timescale is much smaller than the dynami- 
cal timescale since we can take the numerical time step 
being tab < At < tdyn- The computational time of ex- 
plicit and implicit scheme with an equal number of time 
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TABLE 1 
List of Simulation Runs 



model 


K 


r 


state 


P 


Pg 




uv 




K 


RHDSTl 


0.4 


5 
3 


L 


1.0 


3.0 X 10~5 


0.015 








1.0 X 10"* 








R 


2.4 


1.61 X 10""' 


6.25 X 10^3 








2.50 X lO"'^ 


RHDST2 


0.2 


5 
3 


L 


1.0 


4.0 X 10"^ 


0.25 








2.0 X 10-5 








R 


3.11 


0.04512 


0.0804 








3.46 X 10-3 


RHDST3 


0.3 


2 


L 


1.0 


60 


10 








2 








R 


8 


2.34 X 10^ 


1.25 








1.13 X 10^ 


RHDST4 


0.08 


5 
3 


L 


1.0 


6.0 X 10-3 


0.69 








0.18 








R 


3.65 


3.59 X 10-2 


0.189 








1.30 



Note. — Parameter sets of numerical tests. Scattering coefficient cts is taken to be zero in all 
models. 



4.1.2. Radiation transport 

Thi s test has been performed by iTurner fc Stond 
(|2001f) . When the radiation is injected into uniform 
matter, it propagates with the characteristic wave ve- 
locity, while it exchanges energies with the plasma. As 
a test of this effect, we assume that the fluid is static 
and uniform in a simulation box bounded by x — [0, L], 
where L — 1 cm. The radiation is injected from the 
boundary at x = 0. In the uniform medium, the opac- 
ity is set to K = 0.04 cm"^ g. The radiation energy is 
-Elte = 10^° erg cm-'^ and the thermal energy is de- 
termined from the LTE condition. We examined two 
models of the mass density p = 0.25 and 25 g cm-^. 
The corresponding optical depth is r = pnL = 0.01 and 
1, respectively. At the boundary a; = 0, the radiation 
is injected with the energy density of Er = 10^"£'lte- 
The free boundary condition is applied dX x ~ L. Other 
parameters are V = 5/3, and p, = l.O. The Courant- 
Friedrichs-Lewy (CFL) number is taken to be 0.5. 

Figure [5] shows results with r = 0.01. Thick solid 
curves denote the radiation energy density, and thin solid 
curves show the radiation flux at t = 0.5, 1.0, 1.5tc from 
left to right, where Tc = L/c. Since the radiation en- 
ergy is injected from x — the wave front propagates 
from left to right as time goes on. The dashed curves de- 
note the position of the light head Ic = ct/V^. Since we 
assume the Eddington approximation on the closure rela- 
tion, the wave front propagates with the velocity c/VS. 
We can see that the wave front is sharply captured in 
our simulation code. In FLD approximations, the ra- 
diation field evolves obeying the diffusion equation, so 
th at the wave front has a smooth profile (see. Fig. 7 
in iTurner fc: Stonel[2001|) . In our code, although we ap- 
ply the Eddington approximation, the 1st order moment 
equation is solved. Then equations of the radiation field 
have hyperbolic form, so that the wave front can be cap- 
tured using the HLL scheme. 

Figure [3] shows results of the model of larger density 
(r = 1). The radiation propagates with the velocity 

c/\/3. Behind the wave front, the radiation field be- 
comes steady and its energy exponentially decreases with 
x (equivalently, t = npx) due to the absorption. When 
we assume the steady state and the radiation energy is 
much larger than that in LTE (i.e., Er ^ AttB/c), equa- 
tion ([3]) can be solved and we obtain 

cE 

Er — E'o exp(— v^pra), Fr = — ^ exp(— \/3/ckx), (57) 



(|Mihalas fc MihalaslHOSl . Dotted curves in Fig. [2] and 
[3] show solutions obtained from equation ([57)) . We can 
see that numerical results excellently recover analytical 
ones. 

4.1.3. Radiation transport in 2-dimension 

This test has been shown in ITurner fc Stond (|2001[ ) 
that the radiation field propagates in optically thin media 
in X — y plane. The simulation domain is a; = [—L,L] 
and y = [— L, L], where L = 1cm. Numerical grid points 
are {Nx,Ny) = (400,400). The mass density of uniform 
fiuid is p — 0.25 g cm"'^. The other quantities are the 
same as those in § 14.1.21 At the initial state, a larger 
radiation energy is given inside r = y^x'^ + y"^ < 0.1 cm 
and its energy is Er = IO^'^-Elte- The free (Neumann) 
boundary condition is applied on each side of the box. 
The CFL condition is taken as 0.5. 

Figure m shows contours of radiation energy density Er 
at t = 1.5tc. The solid, dashed, and dotted curves de- 
note contours at 10%, 0.1%, 10-'^% of its maximum. Due 
to the enhancement of radiation energy around the ori- 
gin at the initial state, the radiation field propagates in 
a circle, forming a caldera-shaped profile of the radia- 
tion energy density. The wave front propagates with the 
speed of c/-\/3- Figure [5] shows one dimensional cut of 
results at t = 1.5tc. Asterisks, squares and filled circles 
show results along y — 0, x — 0, and y = x, respec- 
tively. We can see that the wave packet has a sharp 
structure. The radiation energy is mainly accumulated 
around the wave front, while it is very small inside the 
wave front (caldera floor). Such the caldera structure 
is not succes sfully reproduced whe n we apply FLD ap- 
proximation (|Turner fc Stondl200lD since the FLD is for- 
mulated based on the diffusion approximation. Although 
the wave speed reduces to c/\/3, solving the 1st order mo- 
ment in radiation field (radiation flux) with Eddington 
approximation has an advantage for studying the prop- 
agation of the radiation pulse in optically thin media. 

In the problem of ii l4.1.2l and ij 14. 1.31 the radiation field 
passes the boundary x = L at t ~ V^Tc- When we adopt 
the free (Neumann) boundary conditions, most of the ra- 
diation energy passes the boundary, while some part of 
it is reflected and stays in the simulation box. An am- 
plitude of the reflected wave i?rof is smaller than that 
passing the boundary E'pass, niax[i5ref/£'pass] ~ 0.8% for 
the one-dimensional test and ~ 5% for the two dimen- 
sional test. The waves are mainly reflected at the corner 
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Fig. 6. — Profiles of p, pg, , E'^ and _F/.^ from top to bottom at 
t = 5000 for the model RHDSTl. Dots and solid curves denote for 
numerical and semi-analytical solutions. 

of the simulation box in the two dimensional simulation 
(i.e., around [x,y] — [±L,±L]). The simple free bound- 
ary condition can be applied in the radiation field with 
this accuracy. 

4.2. Numerical Test for Relativistic Radiation 
Hydrodynamics 

Next, we consider the coupling between the radiation 
and matter fields. We solve RRHD equations with a 
second order accurate scheme. 

Test problems for ra diative shock shown in this subsec- 
tion are developed bv lFarris et al.l ()2008() who proposed 
and solved four shock tube test problems. Initial condi- 
tions of these problems are listed in Table [TJ The initial 
state has a discontinuity at a; = and the system is 
in LTE. Such an initial discontinuity breaks up generat- 
ing waves. When the waves pass away from the simula- 
tion box, the system approaches the steady state. Then 
we compare our numerical results with semi-analytical 
solutions af t er the system reaches a steady state. In 
IFarris et al.l (|2008D . the initial condition is constructed 
by boosting semi-analytical solutions, so that the shock 
moves with an appropriate velocity. In our test, the 
shock is assumed to rest around a; = (shock rest frame). 
Such the problem can be more stringent tests for our code 
to maintain the stationarity (Zanotti ct al. 2011,) . 

The simulation box is bounded by x = [— where 
L = 20 in the normalized unit, and a number of nu- 
merical grid points is Nr^ = 3200. The free boundary 
condition is applied at x = ±L. We give the physi- 
cal quantities at the left and right st ate of p, jp, u, Ej. an d 
flux is assumed to be zero. Following IFarris et all (|2008D . 
the Stcfan-Boltzmann constant is normalized such that 
AiruR = E^j^/Tl = E^ piJTli, where the subscripts L 
and R denote the quantities at the left [x < 0) and right 
{x > 0). 
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Fig. 7. — Profiles of p, pg, , E'^ and F^F from top to bottom at 
t = 5000 for the model RHDST2. Dots and solid curves denote for 
numerical and semi-analytical solutions. 

We note that although we adopt an iteration method 
to integrate stiff source term (equation (|34|)). solutions 
in the following tests converge within the relative error 
of 10~^ without iterations. 

4.2.1. Non-relativistic shock 

In this test, the non-relativistic strong shock exists at 
X = (model RHDSTl). The ratio of the radiation to 
thermal energy in the upstream {x < 0) is 2.2 x 10^"*. We 
take the CFL condition to be 0.9. Fig. [6] shows profile 
at t = 5000. The physical quantities shown in this figure 
are p, pg, , E'^., F^^ from top to bottom. (We note that 
F^^ , which is the radiation flux measur ed in the comoving 
frame, is different from F^ defined in IFarris et al.|[2(308l 
by factor 7). Dots and curves denote for the numerical 
and semi-analytical solutions. 

Since we give an initial condition by a step function at 
X = 0, waves arisen at the discontinuity propagate in the 
ibx— direction (mainly in the -|-x-direction since the up- 
stream is supersonic). After the waves pass the boundary 
at x = ±L, the system reaches the steady state. 

Since the radiation energy is negligible, fluid quantities 
have jumps around x = similar to pure hydrodynam- 
ical shock. The radiation field, on the other hand, has 
a smooth profile in which the radiation energy is trans- 
ported in front of the shock. The radiation energy and 
flux have no discontinuities. The simulation results are 
well consistent with analytical solutions. 

4.2.2. Mildly relativistic shock 

In this test, the mildly-relativistic strong shock exists 
at X = (model RIIDST2). The ratio of the radiation 
to thermal energy in the upstream is 3.3 x 10~^. We 
take the CFL condition to be 0.9. Fig. |7| shows profiles 
aX t = 5000. The physical quantities shown in this fig- 
ure are p, Pg, v^, E'^.F^^ from top to bottom. We can 
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Fig. 8. — Profiles of p, pg, E'^ and F/,^ from top to bottom at 
t = 5000 for the model RHDST3. Dots and solid curves denote for 
numerical and semi-analytical solutions. 
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see that the radiation energy density is no more contin- 
uous but jumps at a; = 0. In the non-relativistic hmit, 
the radiation energy density and its flux are continuous, 
but it is not the case in general. They are no more the 
conserved v ariables at the shock in the relativistic flow 
(| Farris et al. 2008). The solutions obtained in our nu- 
merical simulation (dots) are qualitatively and quanti- 
tatively consistent with semi-analytical solutions (solid 
curves). 

4.2.3. Relativistic shock 

In this test, the highly-relativistic strong shock exists 
at a; ~ (model RHDST3) . The upstream Lorentz factor 
is ~ 10 and the ratio of the radiation to thermal energy is 
3.3 X 10-2. We take the CFL condition to be 0.9. Fig.|8] 
shows profile aX t — 5000. The physical quantities shown 
in this figure are p, Pg, , E'^,Fl!^ from top to bottom. 

After waves generated at the jump {x = 0) passes the 
boundary, the system approaches the steady state. In 
this case, all the physical quantities are continuous. 

The solutions obtained in our numerical simulation ex- 
cellently recover analytical solutions even when the fiow 
speed is highly relativistic. To validate our numerical 
code, we compute the L-\ norms from 

L^[f]^AxY,NAf^-k{x^)l (58) 

i=l 

where fi is the physical quantity at the z-th grid and /a 
is the semi-analytic solution. 
We perform simulations with number of grid points 
= 400,800,1600,3200. In these tests, we use the 
semi-analytic solution as the initial condition. Figure |9] 
shows the L-\ norms of errors in p,u^ ,pg, E'^., and F'.^. 
We can see that all errors converges at second order in 
Ax. 

4.2.4. Radiation dominated shock 

In this test, the mildly-relativistic, radiation domi- 
nated shock exists at a; ~ (model RIIDST4). The ratio 
of the radiation to thermal energy in the upstream is 20. 
We take the CFL condition to be 0.3. Fig. [TOl shows pro- 
file at t = 5000. The physical quantities shown in this 
figure are p, Pg, , E'^.F^.^ from top to bottom. 

After waves generated at the jump [x = 0) pass the 
boundary, the system approaches the steady state. In 
this case, all the physical quantities are continuous and 
their profiles are very smooth since a precursor generated 
from the shock strongly affects the plasma. 

We again find very good agreement with semi- 
analytical solutions even when the radiation energy dom- 
inates the thermal energy. Figure lTTl shows the L-1 norms 
of errors in p, ,pg^ E'^, and -F"^^. In this test, we adopt 
the analytical solu tion as an initial condition following 
iFarris et al.l (|2008[) . We can see that all errors converges 
at second order in Ax. 

5. DISCUSSION 

In the current study, we construct the relativistic hy- 
drodynamic simulation code including the radiation field. 
The magnetic field, which would play a crucial role in 
the relativistic phenomena, is neglected for simplicity. It 
would be quiet simple to include the magnetic field using 
a well-developed numerical scheme. When we extend our 



radiation hydrodynamic code to the radiation magneto- 
hydrodynamic one, the energy momentum tensor is 
replaced by that of the magnetofluids, 

^MHD = {P^ + b>''n'' - V'\r + [pg + 77^^ (59) 

where \f = \u ■ B, [B + {u ■ B)u] /7} is the covariant 
form of the magnetic fields. The magnetic field evolves 
according to the induction equation for the ideal MHD, 

^^{u''b^' -u^'b'') ^0. (60) 

Then the energy momentum conservation equation and 
the induction equation are explicitly integrated using 
the HLL o r higher order appro xi mate Riemann solvers 
(Mignone fc Bodol l2005l [200l IHonkkila fc JanhunenI 
12007; ,Mignone et al.ll2009l ). The source term describing 
the interaction between the matter and the radiation, is 
integrated by applying our proposed scheme. 

Another simplification made in this paper is adopting 
the Eddington approximation that the radiation field is 
assumed to be isotropic in the comoving frame. The 
wave velocity of radiation field has a propagation speed 
of c/y/3- This would become an issue when the fiow 
speed is relativistic (w ~ c). The flow speed, in principle, 
exceeds the phase speed of light mode so that the ra- 
diation energy might accumulate in front of the plasma 
flow. Also the problem appears when we consider the 
magnetofluids. The fast magnetosonic wave speed can 
exceeds the reduced light speed c/-\/3. Another problem 
is that the radiation does not propagate in a straight 
line since we assume that the radiation field is isotropic. 
When the optically thick medium with a finite volume 
is irradiated from one side, there should be a shadow on 
the other side (H aves fc Norman 2003) . Such a shadow 
is no longer formed when we utili ze the Eddington ap- 
proximation (jGonzalez et al.l [20071) . To overcome these 
problems, we should ad mit an anisotropy of radiation 
fiux in a comoving frame (|Kershawlll9"76l : lMinerbollT97l 
lLevermordll984[ ). In our formulation, we do not specify 
the closure relation until S 13.31 Thus we can employ M-1 
closure by replacing matrix A{u) without the other mod- 
ification. The explicit-implicit scheme of the relativistic 
radiation magnetohydrodynamics with the M-1 closure 
will be reported in the near future. 

6. SUMMARY 

We have developed a numerical scheme for solving spe- 
cial relativistic radiation hydrodynamics, which ensures 
the conservation of total energy and flux. The hyperbolic 
term is explicitly solved using an approximate Riemann 
solver, while the source term describing the interaction 
between the matter and the radiation, is implicitly in- 
tegrated using an iteration method. The advantage of 
the implicit scheme is that we can take the numerical 
time step larger than the absorption and scattering time 
scale. This allows us to study the long term evolution of 
the system (typically, the dynamical timescale). 

When integrating the source term, we need to invert 
matrices C and A in our proposed scheme. We have to 
note that the rank of these matrices used in the implicit 
scheme is very small (4x4 for C and 6x6 for A). This 
is because the interaction between the radiation and the 
gas is described by a local nature (source term) when 
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we solve 0th and 1st moments for radiation. This means 
that the numerical code can be easily parallelized and a 
high parallelization efficiency is expected. 

We note that since the matrix C is only 4x4 ma- 
trix, we can invert it analytically. For the matrix A, it 
is relatively small matrix but it is difficult to invert it 
analytically. Thus we decide to use LU decomposition 
method. Since the LU decomposition method can invert 
the matrix directly without iterations, we obtain 
stably. 

We find that the wave front of radiation field can be 
sharply captured using the HLL scheme even when we 
adopt the simple closure relation P^*-' = ^^■'/S (i.e., the 
Eddington approximation). When we adopt the FLD 
approximations, such a sharp wave front cannot be cap- 
tured due to the diffusion approximations. Thus, solving 
the 1st order moment of radiation has an advantage when 
we consider the optically thin medium although the wave 
speed reduces to c/a/3 and the radiation field is isotropic. 

We adopted the iteration scheme to integrate stiff 
source terms. If we need many counts in iteration 
scheme, it causes a load imbalance between each core 
when the code is parallelized. But we have to stress 
that solutions converged only within two steps even if 
At = lO^tab for the problem of radiative heating and 
cooling appeared in § 14.1.11 For the other tests, solution 
converges without iteration. Thus we can expect that the 
load imbalance due to the iteration scheme is not severe 
even in the more realistic problems.^ 

In the previo us papers (jFarris et al.l l200l 
iZanotti et al.l 120 lit ), radiation fields are defined in 
different frames for primitive and conservative variables 
(mixed frame). In such a method, the radiation moment 
equations are simply described. However, this method is 
not suitable for implicit integration. This is because that 
a Lorentz transformation between two frames is needed. 
By the transformation, expressions of the radiation 
fields become quite complex even if the Eddington 
approximation is adopted. Moreover, the radiation 



stress tensor, Pr, is generally non-liner function of the 
radiation energy density and the radiation flux through 
the Eddington tensor. Thus, extension of explicit 
method to implicit method is not straightforward for 
the relativistic radiation hydrodynamics. 

In our method, we treat radiation fields in the labo- 
ratory frame, only. Such a treatment simplifies implicit 
integration. Although we need to invert 6x6 matrix and 
4x4 matrix, our new method is easier than the method 
with mixed frame. Two matrices of C (which is needed 
for implicit integration) and A (which is needed to com- 
pute Pr) are directly inverted using analytical solution 
and LU-decomposition method. Since we need not to 
use some iteration methods to invert them, our method 
is stable and comparatively easy. Our method would be 
quite useful for numerical simulations of relativistic as- 
trophysical phenomena, e.g., black hole accretion disks, 
relativistic jets, gamma-ray bursts, end so on, since both 
the high density region and high velocity region exists, 
and since the radiation processes are play important role 
for dynamics as well as energetics. 
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APPENDIX 
CALCULATING P^' 



In this appendix, we show how to obtain the radiation stress in the observer frame from Er and by assuming the 
Eddington approximation. P^^ is the 3x3 symmetric matrix and there are six unknowns. Since equation (jSip is a 
linear with Er and F^, we need to solve 6th order linear equations, 

Ap = r, (Al) 

where A = {a,,}, p'^ = (P^i, Pf, Pf, P,i2, P,i3, Pf), and = (P", p22^ ^^33^ ^12^ ^13^ ^23) ^he explicit form of 
an is written as 



aii==l + {ui){ui) 
ai2 = {u2){u2) 

ai4 = {ui){u2) 
ai5 = (ui)(M3) 
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ai6 = (W2)(W3) 
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= — ^^-^2 — + -TTT' ^ ^ 

«54 = + TT^' ^^^^^ 



(1+7)^ 1 + 7 



^55 = J- H , _ ' 1 I ' (Ac}^ j 



(A34) 



(M^)(M-0(m)(m) 



062 = 7— — 72 — , (A36) 

(1+7)"^ 1+7 

«63 = — (1^:^)2 — + ^TT' ^ ^ 

«64 = — ^YTW' — + ^T^' ^^^^ 

«65 = — ^YTW' — ^^TV' ^ ^ 

«66 = i + — p^:^)2 — + • (A40) 



In our numerical code, the matrix A is inverted using LJ7-decomposition method. 
Also {dP^^)/{dEr) and {dP^^)/{dF^) are obtained by solving the following equations 



d^_dr^ d^_dr_ 

dEr dEr' dFk dPk' ^^^^^ 



where, 

dB}i 7^ 



(A42) 



dE, 



a^fe --37-^'^"' + w'^i + u'Sl + ^u^u^Uk. (A43) 

Since the matrix A is a function of u and independent of Er and F^, matrix is computed once before calculating 
P'^i,{dPp)/{dEr) and {dP^^)/{dF^). 
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